Type-dependent irreversible stochastic spin models for biochemical reaction networks 
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We describe an approach to model biochemical reaction networks at the level of promotion-inhibition circuitry 
through a class of stochastic spin models that depart from the usual chemical kinetics setup and includes spatial 
and temporal density fluctuations in a most natural way. A particular but otherwise generally applicable choice 
for the microscopic transition rates of the models also makes them of independent interest. To illustrate the 
formalism, we investigate some stationary state properties of the repressilator, a synthetic three-gene network 
of transcriptional regulators that possesses a rich dynamical behaviour. 
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I. INTRODUCTION 

The modeling of biochemical reaction networks is tradi- 
tionally carried out through rate equations based on tech- 
niques inherited from the field of chemical kinetics, some- 
times with refinements such as the use of time-delayed terms, 
differential-difference equations, and stochastic perturbations 
[1-3]. However, the central paradigms of chemical kinet- 
ics, namely, the law of mass-action and the well-stirred re- 
actor approximation, are valid only for slow processes oc- 
curring in dilute solutions at local equilibrium that hardly 
hold in the crowded cellular environment — in a typical cell, 
macromolecules can occupy as much as 40% of the total 
volume available in concentrations of 50^400 mg/ml, with 
steric repulsion effects contributing to the toughness of the 
medium [4]. Chemical master equations, a mesoscopic ap- 
proach to chemical kinetics that possesses connections with 
several branches of nonequilibrium statistical mechanics, are 
also based on the well-stirred approximation [5, 6]. In order 
to circumvent the limitations of the rate equations and also 
to provide modeling tools at varied levels of abstraction, ap- 
proaches based on Boolean networks, stochastic Petri nets, 
and rule-based formalisms, among others, have been devel- 
oped [7]. While some of these approaches are innovative 
in proposing new forms of representing biochemical reaction 
networks and integrating the models with laboratory tools and 
automation (and, as such, sometimes are more of a metamod- 
eling nature), most seldomly abandon chemical kinetics ideas 
for quantitative predictions. 

Here we present an approach to the modeling of biochem- 
ical kinetic phenomena akin to stochastic spin models that 
seems promising [8, 9]. It can in principle represent essential 
features of the components of the systems more directly, pro- 
viding constraints on parameters associated with behaviours 
that can be observed in the wet laboratory, potentially allevi- 
ating the paramenter inference step that greatly hampers semi- 
phenomenologial approaches based on rate equations [10]. 
The spin-like models presented here can also be explored to 
address the practical and difficult question of putting together 
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deterministic kinetics associated with continuous variables 
and stochastic kinetics associated with discrete variables, both 
of which occur in processes of biochemical interest, thus pro- 
viding an alternative to the analysis of biochemical pathways 
where stochasticity is known to play a role [11-16]. 



II. TYPE-DEPENDENT STOCHASTIC SPIN MODELS 

Let T — {di, a.2, . . . , a n } be a finite set of n types (e.g., 

genes or proteins), S a — {s^, s a 2 \ . . . , si S "' ) } the set of S a 
possible internal states of type a, and £ = {(a, s) : a G T, 
s G S a }- Also, let 1/ be the vertex set of a simple graph 
of order V = We call the ordered pair (i, a) G X — 
1/ x 1 a site and denote its internal state by rjf G S a , the 
state space of configurations r/ = (rjf) being given by ft = 
•Sai x So* x ' ' ■ x $a ■ Sites interact through a set of two-body 
interaction matrices Jy ( • , • ) : £ X £ — > R, one for each pair 
of positions i,j G 1A Interactions between types do not need 
to be symmetric, ^ ; otherwise, we will only consider 
isotropic interactions, = JjP, The element J§ , rfy 
denotes the interaction strenght that site (i, a) in the internal 
state r)f exerts upon site (j,b) in the internal state rjj. From 
these matrices we define an "energy" function H : ft — > K as 

(1) 

where Xj is a neighbourhood of (j, b) that may or may not in- 
clude j, b, or (j, b). If rjf promotes r/j, Jij(r)f, rq) < 0, while 
if Til inhibits J?/ : (V? , Vj) > °- Viewed as a spin Hamil- 
tonian, H(rf) is closely related with TV-colour Ashkin-Teller 
and Potts models [17, 18], but generalises them on the counts 
that it is in general a mixed-spins model, since the internal 
state spaces S a do not need to be all identical, and that inter- 
actions between different types do not need to be symmetric. 

Function H (rj) allows us to define a dynamics for the tran- 
sitions of the internal states of the sites from the change in 
H(r]) brought by the transitions, as with the usual stochas- 
tic spin models [8]. Here we will consider single-site tran- 
sitions, although stirring can be added with some extra care. 
Let rjf(s) G f2 be the configuration given by [rjf (s)]j = s if 
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(j,b) = (i,a)md [r 1 f(s)} b 1 



i , — 77* otherwise. The energy cost 
of a transition 77" (r) — > 77? (s) is then given by 



(2) 



Because of the asymmetry in the interactions, A"(r, s)(r]) de- 
composes into Af(r, s)(tj — >■ i) + A"(r, s)(rj <— i), where 



A^(r, S )(r,^z)= ^ [4f(^, S )-J^(^r) 



(3) 



collects the energy difference due to the action of the sites in 
77 upon the site (i, a) when it flips from r/f = r to 77" = s, and 



A?(r,*)fa<-i) = J] 



(4) 



collects the energy diference due to the action of the site (i, a) 
upon the sites of 77 when it flips from 77? — r to r/f = s. 
We now define a dynamics for the model specified by H(rf) 
through the set of single-site transitions rates 



c?(r, «)(!,) = 0(A?(r, *)(!!-►»)). 



(5) 



where 6 : K — > R + is any non-increasing function obeying 
the condition 6(A)e A = e(-A)e~ A . 

The transition rates (5) depend only on the energy differ- 
ence of the single site that flips, not on the global energy dif- 
ference caused by the flip. From the vantage point of the flip- 
ping site, it is as if the rest of the system acted as a heat reser- 
voir that goes unperturbed by the flip — only subsequent flips 
will eventually notice the change. This diverts from the usual 
recipe and has the important consequence that the stationary 
states of the model will not be distributed according to the 
Gibbs measure hg{v) ex P(~H(r])), although there may 
be some function of 77, different from H(r/), that renders a 
Gibbs-like stationary distribution for the model. For reversible 
stochastic spin models, single-site transition rates given by 
cf(r,s)(ri) = 0(A°(r, s)(t?)) guarantee that the stationary 
state will be distributed according to fj,a(r]). For symmetric 
interactions, J° b = jf", we obtain from eqs. (3) and (4) that 
A"(r, s)(tj) = 2A"(r, s)(r) — > i), and the two prescriptions 
coincide up to a factor of 2. So, why should one pick the 
transition rates given by eq. (5) instead of those that guarantee 
that the system will relax to its equilibrium Gibbs distribution? 
The answer is that the rates in eq. (5) lead to forward Kol- 
mogorov equations that, in the mean field approximation — 
corresponding to a well-stirred solution — and in the limit of a 
large number of particles are equivalent to a dynamical system 
x t = V(x t ) for the density profile x t £ K E , with a smooth 



drift vector field V(x t ) 



of the form f(x t )—g(x t 



The rates given by eq. (5) are thus the ones that correctly es- 
tablish the connection between the microscopic description in 
terms of the Markov jump process governed by H (77) and the 
macroscopic description in terms of rate equations. This re- 
sult was obtained in [9] and is mildly related with results first 
obtained by T. Kurtz in the 1970s [19], but the introduction 
of the type-dependent stochastic spin models (1) and the rates 
(5) is novel and provides a versatile modeling framework of 
independent interest. 
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FIG. 1. The repressilator genetic regulatory network circuit. Blunt 
arrows indicate inhibition through a genetic regulation mechanism 
described in the text. 



The simplest type-dependent stochastic spin model has all 
internal state spaces S a = { — 1, +1} and will be referred to 
as type-dependent stochastic Ising model (TDSIM). The most 
general interaction Hj(rj) for TDSIMs is, to within an irrele- 
vant additive constant, given by 



(6) 



where now J°f 



Aab 



and Bf b are scalar quantities. We re- 
mark that Ising-like Hamiltonians have already been used to 
model gene-gene interacting networks, but within the context 
of equilibrium distributions [20]. In our dynamic approach, 
the rates (5) are as important as H b (r]) itself. Notice also 
that the present approach is only barely related with the use 
of Ising spins to analyse consistency and monotonicity of 
reaction network graphs [21], although the determination of 
H b Ari) depends on such graphs. 



III. THE TDSIM FOR THE REPRESSILATOR 

Let us illustrate the formalism by considering the repressi- 
lator, a genetic regulatory network designed to exhibit stable 
oscillations that are believed to be important in the determi- 
nation of the circadian rythms observed in most living organ- 
isms. The repressilator was induced in the prokaryote bac- 
teria Escherichia coli through a genetically engineered plas- 
mid, together with a reporter plasmid that expresses the green 
fluorescent protein (GFP). In this system, the protein LacI 
from E. coli inhibits the transcription of a second gene, tetR 
from the tetracycline-resistance transposon TnlO, whose pro- 
tein product TetR inhibits the transcription of a third gene, cl 
from the A-phage, whose protein CI inhibits the expression 
of lad, closing the loop of negative feedback [22]. This ge- 
netic regulatory network is represented in Figure 1. This is 
clearly a highly stylised description of the true biochemical re- 
action network, that involves different operator sites, depends 
on how many proteins bind to the sites, and have lots of inter- 
mediate steps. It can, however, capture the essential nature of 
the interactions and is widely used to represent biochemical 
networks at a higher level of abstraction. 

The TDSIM for the repressilator in the absence of external 
driving (Afj = Bf b = 0) has three coupling constants, one 
for each pair of unidirectionally interacting types, all positive 
and that can be taken homogeneous. Here we will take all 
coupling constants equal, J AB = J BC = J CA = J, that 
despite being a considerable simplification of the full H b (rf) 
possesses a rich dynamical behaviour already in the mean field 
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approximation [9]. In this case, the two-body interaction term 
becomes 

H (V) = JY1 btrif + vfvf + ViVf] ■ (7) 
The main quantities of interest are the empirical densities 

where <5( • , • ) is the Kronecker delta symbol. For TDSIMs 
we can measure p a — {]~/V)J2i£v 1 li instead, from which 
Pa = |(1 ^ Pa) can be easily recovered. The time evolu- 
tion of these quantities in the stationary state of the model 
for some choices of J appears in Figure 2. All data were 
obtained by Monte Carlo simulations using a heat bath pre- 
scription 9(A) = 1/(1 + e 2A ) for the rates (5) in a simple 
square lattice of V = 100 x 100 sites with periodic bound- 
ary conditions and nearest-neighbour interactions. Notice that 
we include a given position in its own neighbourhood to allow 
for intrasite interactions between different types. One Monte 
Carlo step equals nV move attempts at randomly chosen sites 
(i, a), where n is the number of different types in the system. 

Figure 2 displays the density profiles in the nonequilibrium 
stationary state of the model. From that figure we clearly see 
that the densities of different types oscillate and are out of 
phase. Notice that the curves are mostly pairwise anticorre- 
lated and that different types alternate in the peaks. The os- 
cillations in fig. 2 are similar to the oscillations found experi- 
mentally as well as in ODE models and stochastic simulations 
[22, 23]. When J w 0, the types become independent or 
nearly independent and their densities fluctuate at will, so that 
we do not observe true oscillations. We could identify oscil- 
lations in our finite system for J > 0.07. There is nothing 
special about this value, only that we can clearly observe os- 
cillatory behaviour above it. We found that the amplitudes of 
the oscillations vary little in the range 0.07 < J < 0.42, but 
decay for J > 0.42 and gets smaller as J gets larger past this 
point. We also found that the amplitudes of the oscillations 
scale like yV, signaling that the oscillations are spatially un- 
synchronised, since otherwise the amplitudes would scale like 
V. As a consequence, it becomes difficult to distinguish cy- 
cles or quasi-cycles out of the noise directly from the den- 
sity profiles, and the analysis of correlation functions becomes 
preferable. This is well known from the study of population 
dynamics [24, 25]. We then compute the density-density time 
correlation functions in the stationary state, 

C ab (t) = lim i / [p a (t + t')-p a ][p b (t')-p b }dt', (9) 

and their power spectral densities 

S ab H = — / C ab (t)e- lult dt, (10) 

where p a and p b are the average densities of types a and b in 
the stationary state. In practice, the integration limits in (9) 
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FIG. 2. Evolution of the densities of the types in the stationary state 
of the TDSIM for the repressilator with J — 0.3 (top panel) and 
J — 0.5 (bottom panel). The densities clearly oscillate out of phase 
and are pairwise anticorrelated most of the time. The oscillation am- 
plitudes at J = 0.3 are typical in the whole range 0.07 < J < 0.42. 



and (10) are bounded by the lengths of the time series avail- 
able. In our simulations we sampled the stationary densities 
every At = ^ MCS for 10 4 MCS. 

Figure 3 displays the autocorrelation function Caa^) at 
J = 0.415 normalised by its value at t = and some as- 
sociated Fourier transforms SAA(t)- The other autocorrela- 
tion functions behave like CUx(i) because of the symmetry 
between the types. We see from fig. 3 the decay of the auto- 
correlation function, typical of stochastic dynamics due to the 
variability of the oscillations, and the peak in Saa{u) around 
lu* = 0.26±0.03 MCS" 1 at J = J*. The oscillation frequen- 
cies do not vary much with J as long as J < J* ; otherwise, 
the oscillations cease almost completely for J > J*. 

In Figure 4 we exhibit snapshots of the sites where r/f — 
Vi = Vi m the stationary state for some values of J. This fig- 
ure depicts a typical transition from a disordered phase to an 
antiferromagnetic-like phase. We clearly see how the dynam- 
ics of the types in the stationary state becomes more and more 
constrained by their repressors in the immediate neighbour- 
hood as J gets larger, hence the smaller amplitudes in the os- 
cillations of the densities. From figs. 2 and 4 we can infer that 
there is a transition from a spatially uncorrelated, oscillating 
density stationary state to an almost frozen, non-oscillating 
density stationary state at J* ~ 0.415. The system does not 
freeze completely because of the frustration induced by the 
intrasite interactions between types and the form of the rates 
(5), that depend only on the single site that flips and its neigh- 
bourhood, not on the state of the entire system. We located J* 
by computing the "staggered densities" in lattices of several 
sizes. In the dynamical mean-field approximation to the same 
model this transition could be identified with a Hopf bifurca- 
tion at J* = 2/cos(7r/3) = 4 [9]. We remark that in either 
case the transition at J* should be understood as a change in 
the regime of the dynamical system, not as a thermodynamic 
phase transition, although for systems described by a function 
like H(r]) the two interpretations conflate largely. 

In the actual repressilator, the densities of proteins per cell 
oscillate with an observed period T b s = 160±40 min [22]. In 
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FIG. 3. Autocorrelation function CAA{t) at J = 0.415 normalised 
by its value at t = (upper panel) and some Fourier transforms 
Saa(oj) for several different values of 0.1 < J < 0.415 (lower 
panel). The curve Saa{oj) for J — 0.415 (bolder line) peaks at 
uj* = 0.26 ±0.03 MCS" 1 . 
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FIG. 4. Correlation between the three types in a square lattice of 
100 x 100 sites. The figure depicts the sites with rjf = r)\ — r\1 
(black dots) in the stationary state when J = 0.3 (left panel), J = 
0.415 (mid panel), and J = 0.5 (right panel). At J = 0.5 we see an 
almost exact splitting into two sublattices. In this state, the remaining 
dynamics, responsible for the residual small amplitude oscillations 
shown in the botton panel of fig. 2, occurs mostly in the interstices 
between the sites with "pinned" r/f = rjl = rjl. 



our simulations, we found that at J = J* = 0.415 the period 
Tsim = 3.9 ± 0.4 MCS. We thus have the approximate equiv- 
alence 1 MCS ~ 41 ± 7 min in the real system. Translation 
of these figures into meaningful quantities like transcription 
rates is a delicate question that we intend to pursue elsewhere. 



depicted in fig. 1 with a model description at the same level 
of abstraction. The models can capture several characteristics 
of the system, are predictive, relatively simple, easily com- 
putable, and verifiable in a phenomenological sense. They can 
also be easily composed to describe interacting subsystems, 

H(v,t) = H(v) + H($ +£E *«W>$' (ID 

<j,b) (;,<*) 

in accordance with modularity principles commended by the 
systems approach to biochemical reaction networks [26]. 

We showed that the TDSIM for the repressilator generates 
density oscillations that reproduce those found experimentally 
and in ODE-based models. To display oscillations is a non- 
trivial task for nonequilibrium stationary states and is only 
possible for TDSIMs because the rates (5) do not obey the de- 
tailed balance condition with respect to its "energy" function 
(1) that determines the dynamics. 

The lattice structure of the spin systems provides a natu- 
ral setting to study the spatiotemporal dynamics of extended 
networks, an aspect of biochemical reaction networks that has 
received increasing attention in the context of coupled gene 
regulatory networks [27-32]. Models can include diffusion 
through a Kawasaki-type exchange dynamics and also ac- 
count for the possibility that types may be absent, not only in- 
active, in a given site, e.g., by taking some S a = {— 1 5 0, +1}. 
This possibility allows the modeling of deterministic and 
stochastic kinetics concurrently by putting on the same model 
types of low density (e.g., plasmid copies or enzymes) de- 
scribed by discrete variables r/f together with types of higher 
density (e.g., peptides or small substrate molecules) described 
by an effective density in a mean-field-like description. 

It may be that some biochemical reaction networks give rise 
to TDSIMs resembling Hamiltonians known from other con- 
texts. For example, the circadian oscillations of the proteins 
KaiA, KaiB, and KaiC in cyanobacteria can be modeled by 
the promotion-inhibition circuit A — > C H B — > A [33- 
35], whose TDSIM is closely related with an Ising version 
of the spin-i ferromagnetic-ferromagnetic-antiferromagnetic 
trimerised Heisenberg chain [36], an important model in the 
study of magnetisation processes in strong fields. On the other 
way around, the dynamics of an activator-repressor clock 
model that displays both toggle switch and oscillatory be- 
haviours [37] may be modeled by a dimerised ferromagnetic- 
antiferromagnetic Ising chain that seems unexplored. 

We finally remark that the formalism presented here readly 
applies to non-biochemical reaction networks as well, provid- 
ing a framework in which spatially distributed transformations 
are dealt with in a most natural way. 
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